// ==================================================================================
// Shared Genomics Project MPI Codebase
// Version 1.0 30/04/2010
//
// (c) 2010 University of Manchester all rights reserved
//
// This file is distributed under the GNU General Public License, Version 2.  
// Please see the file COPYING.txt for more details
// ==================================================================================

/*!
\file
\brief	Main Body of the FORMAT_DATA program
*/

/*! 
\defgroup	format_data FORMAT_DATA - Genomic Data-set Filter
\ingroup    mpi_stats
\brief		FORMAT_DATA - Genomic Data-set Filter
\deprecated replaced by the MPI program \ref pformat_data2
\details
<p>This application provide basic filtering for allelic data and is based on the PLINK filtering functions.
It always applies the standard filters of missingness (SNP or Person) and MAF using default threshold or
settings from command-line parameters.
The Application of the Hardy-Weinberg Filter is dependent upon specifying the <em>hwe</em> parameter.</p>

<p>This program also generates randomisation artifacts for the processor cores on a computer cluster.</p>

<p>This program is single-threaded from %main() but 
it can utilise OpenPM to run the filtering codes in more than 1 sub-threads 
if a host machine is multi-core.</p>

<p>If you are using an enterprise release of VC++, this program may not compile as
OpenMP is supported in that release of Visual Studio. To build this program as a standard single-core executable,
simply comment out the HAS_OMP definition Macro.</p>

<h2>OPTIONS</H2>

<table class=options border=0>
<tr><th class=param_header colspan=1>PARAMETER</th>
<th class=param_header colspan=3>DESCRIPTION</th></tr>

<tr><td class=param_name colspan=1>--id STRING</td>
<td class=param_desc colspan=3>Unique identifer for a cluster job (also sets output file stem).</td></tr>

<tr><td class=param_name colspan=1>--rdir PATH</td>
<td class=param_desc colspan=3>Path to remote working directory on central file system.</td></tr>

<tr><td class=param_name colspan=1>--ldir PATH</td>
<td class=param_desc colspan=3>Path to local working directory on each core of the cluster.</td></tr>

<tr><td class=param_name colspan=1>--bmap FILE</td>
<td class=param_desc colspan=3>Path to the BMAP file (binary-encoded SNP list).</td></tr>

<tr><td class=param_name colspan=1>--bped FILE</td>
<td class=param_desc colspan=3>Path to the BPED file (binary-encoded pedigree data).</td></tr>

<tr><td class=param_name colspan=1>--phe FILE</td>
<td class=param_desc colspan=3>
Path to  text file containing the phenotype data. This file accessed by index or by named columns.
Missing phenotype values labelled as -9.
</td></tr>

<tr><td class=param_name colspan=1>--mpheno 1..n</td>
<td class=param_desc colspan=3>Phenotype column index to read from text file. Indices start from '1' (cf. PLINK)</td></tr>

<tr><td class=param_name colspan=1>--pheno-name STRING</td>
<td class=param_desc colspan=3>A named column to read from the phenotype text file.</td></tr>

<tr><td class=param_name colspan=1>--hwe 0.0..1.0</td>
<td class=param_desc colspan=3>Apply Hardy-Weinberg %filter and specify the threshold</td></tr>

<tr><td class=param_name colspan=1>--mind 0.0..1.0</td>
<td class=param_desc colspan=3>Set the threshold person missingness %filter</td></tr>

<tr><td class=param_name colspan=1>--geno 0.0..1.0</td>
<td class=param_desc colspan=3>Set the threshold SNP missingness %filter</td></tr>

<tr><td class=param_name colspan=1>--maf 0.0..1.0</td>
<td class=param_desc colspan=3>Set the threshold Minor allelic frequency (MAF) %filter </td></tr>

<tr><td class=param_name colspan=1>--available-procs 1..n</td>
<td class=param_desc colspan=3>Set number of processor cores in the computer cluster</td></tr>

<tr><td class=param_name colspan=1>--row-major i or s </td>
<td class=param_desc colspan=3>Set the row-major of the filtered data-set. 
Either <em>i</em> for Individual-major or <em>s</em> for SNP-major.
</td></tr>

<tr><td class=param_name colspan=1>--nperm 1..10000 </td>
<td class=param_desc colspan=3>Set number of permutations if need to generate permutation index file to control
%sample randomisation.
</td></tr>

<tr><td class=param_name colspan=1>--seed NUMBER</td>
<td class=param_desc colspan=3>Set the seed value for the PRNG. 
If this value not specified at the command-line, defaults to the current system clock.
</td></tr>

<tr><td class=param_name colspan=1>--exclude-list FILE</td>
<td class=param_desc colspan=3>Path to text file containing a person exclusion list, i.e. people
to remove from a data-set. List simply indices of people contained in the %sample array.
Text formatted, 1 number per line.
</td></tr>

<tr><td class=param_name colspan=1>--missingness FILE</td>
<td class=param_desc colspan=3>Path to pre-calculated person missingness file.
</td></tr>

<tr><td class=param_name colspan=1>--founders FILE</td>
<td class=param_desc colspan=3>Path to founders file containing gender data of people in the %sample array.
</td></tr>

<tr><td class=param_name colspan=1>--no-rng</td>
<td class=param_desc colspan=3>If state flag set, suppresses the creation of randomisation artifacts like the seed file.
</td></tr>
</table>
*/

/*! 
\brief Flag to state if compiler supports OMP 
\ingroup format_data
*/

//#define HAS_OMP

#ifdef HAS_OMP
#include <omp.h>
#endif

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#include <math.h>

#include "SMG_All.h"
#include "bitstring.h"

#define FILTER_MODE 1
#define WRITE_MODE 2
#define CANCEL_MODE 3
#define NREQUIRED_ARGS 9

int format_main(int argc, char *argv[]);
int main(int argc, char *argv[]);
void cleanUp();
int snp_do_maf_and_mind(char *bpedData, struct small_locus *snps, int nSNPs, BOOL *snps_to_del, char *snps_reasons_to_del, double *snp_del_stats, int nSamples, BOOL *persons_to_del, BOOL *sex, BOOL *founders, struct chromosomes *chromos, struct selected_options *ops);
int snp_do_hwe(char *bpedData, struct small_locus *snps, int nSNPs, BOOL *snps_to_del, char *snps_reasons_to_del, double *snp_del_stats, int nSamples, BOOL *persons_to_del, BOOL *sex, BOOL *founders, double *phe_values, struct chromosomes *chromos, struct selected_options *ops);

/*!
\brief Command-line options of FORMAT_DATA
*/
struct option long_options[] = {
	{ "bmap", required_argument, NULL, 'm' }, // Path to the BMAP file.
	{ "bped", required_argument, NULL, 'B' }, // Path to BPED File.
	{ "phe", required_argument, NULL, 'f' }, // PHE phenotype file.
	{ "id", required_argument, NULL, 'i' }, // Run ID
	{ "ldir", required_argument, NULL, 'd'}, // Output direct
	{ "rdir", required_argument, NULL, 'R'}, // Output direct
	{ "mpheno", required_argument, NULL, 'n'}, // Index number of data column in the PHE file.
	{ "hwe", required_argument, NULL, 'a'}, // Default Hardy-Weinburg filter
	{ "hwe2", required_argument, NULL, 'b'}, // The Other Hardy-Weinburg filter
	{ "mind", required_argument, NULL, 'e'}, // Missingness per individual.
	{ "geno", required_argument, NULL, 'g'}, // Missingness per marker
	{ "maf", required_argument, NULL, 'h'}, // Minor Allele Frequency
	{ "available-procs", required_argument, NULL, 'N'}, // Number processors in cluster, used to size the random seed file.
	{ "row-major", required_argument, NULL, 'r'}, // Row major of PED file (SNP or individual).
	{ "nperm", required_argument, NULL, 'W'}, // No. of permutations
	{ "seed", required_argument, NULL, 'S' }, // Seed value for the random number generator.
	{ "mode", required_argument, NULL, 'M' }, // mode to run this programme.
	{ "exclude-list", required_argument, NULL, 'E'}, // Sample exclusion list a.k.a category filter.
	{ "missingness", required_argument, NULL, 'Q'}, // Path to sample missingness file.
	{ "founders", required_argument, NULL, 'F'}, // Path to founders/gender data file.
	{ "no-rng", optional_argument, NULL, 'z' },  // Suppress the generation of the RNG files (compatibility flag for solcol.exe).
	{ "pheno-name", required_argument, NULL, 'Z' },  // pheno column name.
	{ 0,    0,    0,    0   } 
};

struct small_locus *map = NULL;
struct chromosomes _chromos;
struct pedfile _p;
struct mapfile _m;
double *_phe_values = NULL, *missingness_per_person = NULL, *filter_statistics_local = NULL, *filter_statistics_global = NULL;
BOOL *founders_data = NULL, *_persons_to_del = NULL, *snps_to_del_local = NULL, *phe_missing = NULL, *snps_to_del_global = NULL;
char *_bpedData = NULL, *reasons_to_del_local = NULL, *reasons_to_del_global = NULL;

/*!
\brief Main body of the FORMAT_DATA program
\ingroup format_data
*/
int format_main(int argc, char *argv[]) {
	char ch;
	struct selected_options ops;
	int nRequiredArgs = NREQUIRED_ARGS, mode = FILTER_MODE, nAvailableProcs = AVAILABLE_CORES;
	BOOL isSnpMajor = FALSE, isIndividualMajor = FALSE, jobDone = FALSE, hasExcludeList = FALSE, do_rng = TRUE;
	int myrank = 0, numprocs = 1;
	char rowmode = 's';
	unsigned long seed = (unsigned long) time(NULL); 
	int nSNPs = 0, nMap = 0, nSamples = 0;

	// Initialise some stuff.
	options_init(&ops);
	pedfile_init(&_p);
	mapfile_init(&_m);
	chromosomes_init(&_chromos);
	chromosomes_define_human(&_chromos);
	*ops.ldir = *ops.rdir = '\0' ; 
	ops.verbose = FALSE;

	if (argc == 1) return -1;

	while ((ch = (char) getopt_long(argc, (const char**) argv, "m:v:p:f:i:", long_options, NULL)) != -1) {
		switch (ch) {
			double value;
			unsigned long _seed;

			case 'm': // BMap File
				strcpy(ops.szBmapFileName, optarg);
				nRequiredArgs--;
				break;
			case 'B': // Binary BPED File
				strcpy(ops.szBpedFileName, optarg);
				nRequiredArgs--;
				break;
			case 'f': // Phenotype file
				strcpy(ops.szPheFileName, optarg);
				nRequiredArgs--;
				break;
			case 'i': // Run ID
				strcpy(ops.szRunId, optarg);
				break;
			case 'd': // Local working directory.
				strcpy(ops.ldir, optarg);
				nRequiredArgs--;
				break;
			case 'R': // Remote working directory.
				strcpy(ops.rdir, optarg);
				nRequiredArgs--;
				break;
			case 'n': // Column to analyse in PHE file
				ops.mpheno = atoi(optarg);
				if (ops.mpheno < 0) ops.mpheno = 1;
				break;
			case 'N': // Theoretical available processors, need to size seed file.
				nAvailableProcs = atoi(optarg);
				if (nAvailableProcs == 0) nAvailableProcs = AVAILABLE_CORES;
				break;
			case 'r': // Row-major mode for the binary formatted PED file.
				rowmode = *optarg;
				if (rowmode == 'i' || rowmode == 's') {
					if (rowmode == 's') { // SNP-major
						isSnpMajor = TRUE; 
						isIndividualMajor = FALSE;
					}  else if (rowmode == 'i') { // individual-major
						isSnpMajor = FALSE; 
						isIndividualMajor = TRUE;
					}
					nRequiredArgs--;
				}
				break;
			case 'a': // Do the HW filter on the population
				ops.HWD_test = TRUE;
				ops.HWD_standard = FALSE;
				value = atof(optarg);
				if (!(value >= 0.0 && value <= 1.0)) value = DEFAULT_HWE_LIMIT;
				ops.HWD_limit = value;
				break;
			case 'b': // Type of HW test and the filter limit.
				ops.HWD_test = TRUE;
				ops.HWD_standard = TRUE;
				value = atof(optarg);
				if (!(value >= 0.0 && value <= 1.0)) value = DEFAULT_HWE_LIMIT;
				ops.HWD_limit = value;
				break;
			case 'e': // Missing limit for Person's SNP missingness.
				value = atof(optarg);
				if (!(value >= 0.0 && value <= 1.0)) value = DEFAULT_MIND;
				ops.MAX_IND_MISSING = value;
				break;
			case 'g': // Limit for SNP missingness
				value = atof(optarg);
				if (!(value >= 0.0 && value <= 1.0)) value = DEFAULT_GENO;
				ops.MAX_GENO_MISSING = value;
				break;
			case 'h': // Major allelic Frequency.
				value = atof(optarg);
				if (!(value >= 0.0 && value <= 1.0)) value = DEFAULT_MAF;
				ops.min_af = value;
				break;
			case 'W':
				ops.nperm = atoi(optarg);
				if (ops.nperm <= 0) ops.nperm = DEFAULT_NPERM;
				break;
			case 'S':
				_seed = atol(optarg);
				if (_seed != 0) seed = _seed;
				break;
			case 'M': // Calling mode for this program (filter or right output).
				mode = atoi(optarg);
				if (mode == FILTER_MODE || mode == WRITE_MODE || mode == CANCEL_MODE) nRequiredArgs--;
				break;
			case 'E': // Sample Exclusion List
				strcpy(ops.szExcludeListFileName, optarg);
				hasExcludeList = TRUE;
				break;
			case 'Q': // path to the 'person' missingness file.
				strcpy(ops.szMissingnessFileName, optarg);
				nRequiredArgs--;
				break;
			case 'F': // path to the founders/gender file.
				strcpy(ops.szFoundersFileName, optarg);
				nRequiredArgs--;
				break;
			case 'z': // Suppress generation of the RNG file (perm index & seed files)
				do_rng = FALSE;
				break;
			case 'Z': // Pheno column-name to read from the PHE file.
				strcpy(ops.szPheColumnName, optarg);
				break;
			default:
				break;
		} 
	}
	if (nRequiredArgs != 0) return -1;

	init_genrand(seed);
	chromosomes_define_human(&_chromos);

	if (mode == FILTER_MODE) {
		if ((_bpedData = bpedfile_load(&ops, &nSNPs, &nSamples)) != NULL) {
			if (_bpedData != NULL) {
				// Allocate the result arrays associated with the filtering.
				_persons_to_del = create_bool_array(nSamples, FALSE);
				phe_missing = create_bool_array(nSamples, FALSE);
				snps_to_del_local = create_bool_array(nSNPs, FALSE);
				reasons_to_del_local = create_char_array(nSNPs, FILTER_REASON_NOTHING);
				filter_statistics_local = create_vector(nSNPs, -1);
				missingness_per_person = missingnessfile_load(&ops, &nSamples);
				founders_data = foundersfile_load(&ops, &nSamples);
				map = bmapfile_load(&ops, &nMap);
				_phe_values = phefile_load(&ops, nSamples, phe_missing);

				if (_bpedData != NULL && _persons_to_del != NULL && snps_to_del_local != NULL && reasons_to_del_local != NULL &&
					filter_statistics_local != NULL && phe_missing != NULL && missingness_per_person != NULL && map != NULL && _phe_values != NULL &&
					nSNPs == nMap) 
				{
					BOOL *sex = founders_data, *founders = founders_data + nSamples; // Set the gender & founder pointers as in the same array.

					if (hasExcludeList) bfilter_parse_exclude_list(&ops, nSamples, _persons_to_del); // Category filter.

					bfilter_person_do_mind(nSamples, missingness_per_person, _persons_to_del, &ops);
					snp_do_maf_and_mind(_bpedData, map, nSNPs, snps_to_del_local, reasons_to_del_local, filter_statistics_local, nSamples, _persons_to_del, sex, founders, &_chromos, &ops);
					if (ops.HWD_test) snp_do_hwe(_bpedData, map, nSNPs, snps_to_del_local, reasons_to_del_local, filter_statistics_local, nSamples, _persons_to_del, sex, founders, _phe_values, &_chromos, &ops);

					// If the BPED write fails there's really no-point carrying on from here.
					if (bpedfile_write_filtered(_bpedData, nSNPs, nSamples, snps_to_del_local, _persons_to_del, myrank, &ops, ops.ldir, ops.rdir)) {
						if (boolfile_create_and_copy(snps_to_del_local, nSNPs, &ops, SUFFIX_SNP_TO_DEL) && boolfile_create_and_copy(_persons_to_del, nSamples, &ops, SUFFIX_PERSON_TO_DEL)) {
							if (pformat_write_filter_log(nSNPs, map, snps_to_del_local, reasons_to_del_local, filter_statistics_local, &ops, ops.ldir, ops.rdir, nSamples, _persons_to_del, missingness_per_person)) {
								if (pformat_write_filter_dims(nSNPs, nSamples, numprocs, snps_to_del_local, _persons_to_del, &ops)) {
									bfilter_write_log(&ops);
									jobDone = TRUE;
								}
							}
						}
					}
				}
			}
		}
	}

	cleanUp();

	if (jobDone) return 0;
	else return -1;
}

int main(int argc, char *argv[]) {
	return format_main(argc, argv);
}

void cleanUp() {
	// Free-up the dynamic memory.
	chromosomes_free_attributes(&_chromos);
	if (map != NULL) { free(map); map = NULL; }
	if (_phe_values != NULL) { free(_phe_values); _phe_values = NULL; }
	if (_bpedData != NULL) { free(_bpedData); _bpedData = NULL; }
	if (missingness_per_person != NULL) { free(missingness_per_person); missingness_per_person = NULL; }
	if (founders_data != NULL) { free(founders_data); founders_data = NULL; }
	if (_persons_to_del != NULL) { free(_persons_to_del); _persons_to_del = NULL; }
	if (snps_to_del_local != NULL) { free(snps_to_del_local); snps_to_del_local = NULL; }
	if (reasons_to_del_local != NULL) { free(reasons_to_del_local); reasons_to_del_local = NULL; }
	if (filter_statistics_local != NULL) { free(filter_statistics_local); filter_statistics_local = NULL; }
	if (phe_missing != NULL) { free(phe_missing); phe_missing = NULL; }
	if (snps_to_del_global != NULL) { free(snps_to_del_global); snps_to_del_global = NULL; }
	if (filter_statistics_global != NULL) { free(filter_statistics_global); filter_statistics_global = NULL; }
	if (reasons_to_del_global != NULL) { free(reasons_to_del_global); reasons_to_del_global = NULL; }
	pedfile_free_attributes(&_p);
	mapfile_free_attributes(&_m);
}

int snp_do_maf_and_mind(char *bpedData, struct small_locus *snps, int nSNPs, BOOL *snps_to_del, char *snps_reasons_to_del, double *snp_del_stats, int nSamples, BOOL *persons_to_del, BOOL *sex, BOOL *founders, struct chromosomes *chromos, struct selected_options *ops) {
	int i, nActualSamples = 0;
	size_t offset = 0;

	if (bpedData == NULL || snps == NULL || nSNPs <= 0 || snps_to_del == NULL || snps_reasons_to_del == NULL || snp_del_stats == NULL ||
		nSamples <= 0 || persons_to_del == NULL || sex == NULL || founders == NULL ||
		chromos == NULL || ops == NULL) return 0;
	else if (!chromos->isdefined) return 0;

	// Count the actual number of samples, i.e. people flagged out to be deleted.
	// If zero, flag all the SNPs for deletion as no more information in the input data-set.
	for (i = 0; i < nSamples; i++) if (!persons_to_del[i]) nActualSamples++;
	if (nActualSamples == 0) {
		for (i = 0; i < nSNPs; i++) {
			snps_to_del[i] = TRUE;
			snps_reasons_to_del[i] = FILTER_REASON_PERSON_MISSINGNESS;
		}
		return 1;
	}

	offset = nSamples * 2 * sizeof(char);
	#ifdef HAS_OMP
	#pragma omp parallel for private(i) shared(nSNPs, bpedData, ops, nActualSamples, nSamples, offset, snps, snps_to_del, chromos, snps_reasons_to_del, snp_del_stats, persons_to_del, sex)
	#endif
	for (i = 0; i < nSNPs; i++) {
		char *pBpedData = bpedData + (i * offset);
		register int j;
		int geno_nm = 0; // count 1 per genotype, for missingness
		int geno_real = ops->oblig_missing ? 0 : nActualSamples; // count 1 per non-obligatory missing genotype (or set to N individuals)
		double snp_genotyping = 0.0;
		BOOL X = FALSE;
		BOOL haploid = FALSE;
		int nm = 0; // count 1 per allele, for frequency
		struct small_locus *loc = &snps[i];
		double bp = 0, pos = 0, freq = 0;

		if (snps_to_del[i]) continue;

		// Determine type of SNP
		if (chromos->chr_sex[loc->chr]) X = TRUE;
		else if (chromos->chr_haploid[loc->chr]) haploid = TRUE;
		
		for (j = 0; j < nSamples; j++, pBpedData += 2) {
			BOOL s1, s2;

			if (persons_to_del[j]) continue;
			s1 = pBpedData[0];
			s2 = pBpedData[1];

			// Check female Y genotypes
			if (chromos->chr_Y[loc->chr] && !sex[j]) {
				// Set to missing, unless in a RECODE mode.
				if (!ops->preserve_all_genotypes) {
					s1 = pBpedData[0] = TRUE;
					s2 = pBpedData[1] = FALSE;
					continue; // Skip to next person.
				}
			}
			
			// For haploid heterozygosity check, also consider all individuals
			if (haploid || (X && sex[j])) {
				if (!s1 && s2) {
					// Set to missing unless in RECODE mode
					if (!ops->preserve_all_genotypes) {
						pBpedData[0] = TRUE;
						pBpedData[1] = FALSE;
					}
				}
			}
			
			// For missing genotypes
			if (!(s1 && !s2)) geno_nm++;

			// For allele frequencies only consider founders.
			if (ops->summ_nonfounders || founders[j]) {
				if (haploid || (X && sex[j])) {
					// Haploid counts
					// "1" allele 
					if (!s1 && !s2) {  // FF = hom(11)
						freq++;
						nm++;
					} else if (s1 && s2) { // TT = hom(22)
						nm++;
					}
				} else {
					// Autosomal count
					// "1" allele 
					if (!s1) { 
						if(!s2) {	// 00 = hom(11)
							freq += 2;
							nm += 2;
						} else {  // 01 = het(12)
							freq += 1;
							nm += 2;
						}
					} else if (s2) {
						nm += 2;
					}
				}
			} 
		}

		// Calculate allele frequencies
		if (ops->af_count) { // Allele counts...
			bp = (long int) (nm - freq);            // Use freq to store count (keep as is)	
			pos = geno_real - geno_nm; // Use "pos" to store number missing genotypes
		} else {
			// ... or frequencies
			if (nm > 0) freq /= nm; 
			else freq = 1;
		}

		// Record total proportion of missingness.
		snp_genotyping =  nActualSamples > 0 ? (double) geno_nm/ (double) geno_real : 0;
		// Exclude if SNP has too many missing genotypes
		if (snp_genotyping < (1- ops->MAX_GENO_MISSING)) {
			snps_to_del[i] = TRUE;
			snps_reasons_to_del[i] = FILTER_REASON_SNP_MISSINGNESS;
			snp_del_stats[i] = 1 - snp_genotyping;
		}

		// Make allele1 always the least common allele
		if (!ops->af_count && freq > 0.5) {
			char tmp[ALLELE_LENGTH];
			char *pBpedRecord = bpedData + (nSamples * 2 * sizeof(char) * i); // Go to the current SNP record & adjust the ordering.

			// then we need to swap alleles
			freq = 1 - freq;
			strcpy(tmp, loc->allele2);
			strcpy(loc->allele2, loc->allele1);
			strcpy(loc->allele1, tmp);

			// Update the data in the BPED array.
			for (j = 0; j < nSamples; j++, pBpedRecord += (2 * sizeof(char))) {
				if (persons_to_del[j]) continue;
				if (pBpedRecord[0] == pBpedRecord[1]) {
					pBpedRecord[0] = !pBpedRecord[0];
					pBpedRecord[1] = !pBpedRecord[1];
				}
			}	
		} 
		
		// Remove rare SNPs
		// Note epsilon correction for MAF, due to floating point 
		// issues: only apply to the lower MAF range
		if (freq < 0 || (freq + ops->epsilon) < ops->min_af || freq > ops->max_af) {
			snps_to_del[i] = TRUE;
			snps_reasons_to_del[i] = FILTER_REASON_SNP_MAF;
			snp_del_stats[i] = freq;
		}
	}

	return 1;
}

int snp_do_hwe(char *bpedData, struct small_locus *snps, int nSNPs, BOOL *snps_to_del, char *snps_reasons_to_del, double *snp_del_stats, int nSamples, BOOL *persons_to_del, BOOL *sex, BOOL *founders, double *phe_values, struct chromosomes *chromos, struct selected_options *ops) {
	int i;
	size_t offset = 0;
	
	if (bpedData == NULL || snps == NULL || nSNPs <= 0 || snps_to_del == NULL || snps_reasons_to_del == NULL || snp_del_stats == NULL ||
		nSamples <= 0 || persons_to_del <= 0 || sex == NULL || founders == NULL ||
		chromos == NULL || ops == NULL || phe_values == NULL) return 0;
	else if (!chromos->isdefined) return 0;

	// Loo through the SNPs
	offset = nSamples * 2 * sizeof(char);
	#ifdef HAS_OMP
	#pragma omp parallel for private(i) shared(nSNPs, bpedData, ops, nSamples, offset, snps, snps_to_del, chromos, snps_reasons_to_del, snp_del_stats, persons_to_del, sex, phe_values, founders)
	#endif
	for (i = 0; i < nSNPs; i++) {
		char *pBpedData = bpedData + (i * offset);
		register int j;
		struct small_locus *loc = &snps[i];
		
		// Compute p-values for HWE test in cases, controls & all
		// Only consider founders
		int a11 = 0, a12 = 0, a22 = 0;
		int u11 = 0, u12 = 0, u22 = 0;
		int b11 = 0, b12 = 0, b22 = 0;

		// Allele frequencies
		double afreq = 0, ufreq = 0, freq = 0;
		BOOL include_cases = TRUE, include_controls = TRUE;
		BOOL X = FALSE, haploid = FALSE;

		// Already removed by some other filter so skip to next SNP record.
		if (snps_to_del[i]) continue; 
		
		if (chromos->chr_sex[loc->chr]) X = TRUE;
		else if (chromos->chr_haploid[loc->chr]) haploid = TRUE;
		
		// Iterate over each individual, founders only
		for (j = 0; j < nSamples; j++, pBpedData += 2) {
			if (persons_to_del[j]) continue;
			
			// Only consider founders, & diploid genotypes
			if (ops->summ_nonfounders || founders[j]) {
				if ( !(haploid || (X && sex[j]))) {	
					BOOL s1 = pBpedData[0], s2 = pBpedData[1];
					
					// Consider everybody, irrespective of phenotype
					// (QT, C/C or missing)
					if (!s1) { 
						if (!s2) b11++;   //   00 = hom(11)
						else b12++;       //   01 = het(12)
					}
					else if ( s2 ) b22++; // 11 = hom(22)

					// for binary trait, separately for cases/controls
					if (ops->bt)  {
						if (phe_values[j] == 1) {
							if (!s1) { 
								if (!s2) u11++;   //   00 = hom(11)
								else u12++;       //   01 = het(12)
							}
							else if (s2) u22++; //   11 = hom(22)

						}
						else if (phe_values[j] == 2) {
							if (!s1) { 
								if (!s2) a11++;   //   00 = hom(11)
								else a12++;         //   01 = het(12)
							}
							else if ( s2 ) a22++; //   11 = hom(22)
						}
					}
				}
			}
		}
		
		// Allele frequencies
		if (ops->qt)
			freq = ( b11 + (double)b12/2.0 ) / (double) (b11 + b12 + b22);
		else {
			afreq = (a11 + (double) a12/2.0) / (double) (a11 + a12 + a22);
			ufreq = (u11 + (double) u12/2.0) / (double) (u11 + u12 + u22);
			freq =  (b11 + (double) b12/2.0) / (double) (b11 + b12 + b22);

			if (a11 + a12 + a22 == 0 ) include_cases = FALSE;
			if (u11 + u12 + u22 == 0 ) include_controls = FALSE;
		}

		if (ops->qt) {
			double p;

			if (ops->HWD_standard) {
				double tot = b11 + b12 + b22;
				double exp_11 = freq * freq * tot;
				double exp_12 = 2 * freq * (1-freq) * tot;
				double exp_22 = (1-freq) * (1-freq) * tot;

				double chisq = ( (b11-exp_11)*(b11-exp_11) ) / exp_11 
					+ ( (b12-exp_12)*(b12-exp_12) ) / exp_12 
					+ ( (b22-exp_22)*(b22-exp_22) ) / exp_22 ;

				p = chiprobP(chisq,1);
			}
			else
				p = SNPHWE( b12, b11, b22 );

			if (p <= ops->HWD_limit ) {
				snps_to_del[i] = TRUE; // Set to delete SNP
				snps_reasons_to_del[i] = FILTER_REASON_SNP_HWE;
				snp_del_stats[i] = p;
			}
		} 
		else
		{
			// For case/control data
			double p, p_a, p_u;

			if (ops->HWD_standard) {
				double exp_a11 = afreq * afreq * (a11+a12+a22);
				double exp_a12 = 2 * afreq * (1-afreq) * (a11+a12+a22);
				double exp_a22 = (1-afreq) * (1-afreq) * (a11+a12+a22);

				double exp_u11 = ufreq * ufreq * (u11+u12+u22);
				double exp_u12 = 2 * ufreq * (1-ufreq) * (u11+u12+u22);
				double exp_u22 = (1-ufreq) * (1-ufreq) * (u11+u12+u22);

				double exp_11 = freq * freq * (b11+b12+b22);
				double exp_12 = 2 * freq * (1-freq) * (b11+b12+b22);
				double exp_22 = (1-freq) * (1-freq) * (b11+b12+b22);

				double chisq_a = ( (a11-exp_a11)*(a11-exp_a11) ) / exp_a11 
					+ ( (a12-exp_a12)*(a12-exp_a12) ) / exp_a12 
					+ ( (a22-exp_a22)*(a22-exp_a22) ) / exp_a22 ;

				double chisq_u = ( (u11-exp_u11)*(u11-exp_u11) ) / exp_u11 
					+ ( (u12-exp_u12)*(u12-exp_u12) ) / exp_u12 
					+ ( (u22-exp_u22)*(u22-exp_u22) ) / exp_u22 ;

				double chisq = ( (b11-exp_11)*(b11-exp_11) ) / exp_11 
					+ ( (b12-exp_12)*(b12-exp_12) ) / exp_12 
					+ ( (b22-exp_22)*(b22-exp_22) ) / exp_22 ;

				p = chiprobP(chisq,1);
				p_a = chiprobP(chisq_a,1);
				p_u = chiprobP(chisq_u,1);
			} else {
				p = SNPHWE( b12, b11, b22 );
				p_a = SNPHWE( a12, a11, a22 );
				p_u = SNPHWE( u12, u11, u22 );
			}

			// Controls (and, if possible, exclude on this value)
			if (include_controls && p_u < ops->HWD_limit && p > -1) {
				if (!ops->HWD_filter_on_all) { 
					snps_to_del[i] = TRUE;
					snps_reasons_to_del[i] = FILTER_REASON_SNP_HWE;
					snp_del_stats[i] = p_u;
				}
			}

			// In total sample, and if needed, exclude here
			if (p < ops->HWD_limit && p > -1) { 
				if (ops->HWD_filter_on_all || !include_controls) {
					snps_to_del[i] = TRUE;
					snps_reasons_to_del[i] = FILTER_REASON_SNP_HWE;
					snp_del_stats[i] = p;
				}
			}
		}
	}

	return 1;
}